22 years of satellite imagery reveal a major destabilization structure at Piton de la Fournaise

Volcanic activity can induce flank failure, sometimes generating large earthquakes and tsunamis. However, the failure structures have never been fully characterized and the failure mechanism is still debated. Magmatic activity is a possible trigger, either through fault slip, which might be induced by dyke intrusions, or through sill intrusions, which might be undergoing coeval normal displacements and slip. At the Piton de la Fournaise volcano, satellite imagery combined with inverse modeling highlights the pathways of 57 magmatic intrusions that took place between 1998 and 2020. We show that a major arcuate dyke intrusion zone is connected at depth to a sill intrusion zone, which becomes a fault zone towards the sea, forming a spoon-shaped structure. Some sills are affected by coeval normal displacement and seaward slip. Overall, the structure is characterized by a continuum of displacement from no slip, to sheared sills and finally pure slip. Repeated intrusions into this spoon-shaped structure could trigger catastrophic collapses.


Inverse modeling methodology Modeling and inversions
In the inverse modeling approach, observation data are used to infer the values of the parameters characterizing a system [7]. Mathematically, the relation between the observations and the model parameters is given by: where d is the vector of observed values, F is the forward modeling operator, which is a function of the input parameters m, and ϵ is the vector of measurements and modeling errors.
For ground deformation, d is the observed displacement and F (m) is a non-linear function of the source location, geometry and pressure. When inverting deformation data derived from Interferometric Synthetic Aperture Radar (InSAR), measurement errors mainly comes from atmospheric effects [8] and discrepancies in the time periods covered by the data, while the modeling errors come from inexact assumption about the source types, number of sources [9,5], inadequate evaluation of the source interactions [10], or an inexact mechanical framework [11,12,13].
Inversion of geodetic data is an ill-posed problem with non-unique solution, leading to different set of parameter values explaining the data equally well. In consequence, any a priori constraints on the models (as the geological significance for instance) are valuable to restrict the range of possible parameter values and better determine a relevant model.
During inversions, the parameter space is explored to determine the range of parameter values that explain the data within their uncertainties by minimizing the difference between the data and the model. It is defined by the cost function: where u obs and u mod are respectively the vector of observed and modeled displacements, and C d is the covariance matrix of the data allowing to weight data according to their variance.
We compute the covariance matrix in order to take into account the correlated random noise of InSAR data [7]. As proposed by Fukushima et al. [2], we assumed an exponential decrease of noise correlation with distance. It is described by: where d is the distance between two points, σ 2 is the noise variance, and a is the correlation length. Variance comes from the combined effects of (1) interferogram uncertainties (estimated in underforming areas) and (2) the uncertainties of the model. Total variance is estimated through the residuals of preliminary inversions. We used a variance value of 5×10 −4 m 2 and a correlation distance of 850 m previously estimated at Piton de la Fournaise [2,6]. Moreover, as ascending and descending interferograms show very contrasted displacement fields at Piton de la Fournaise, it can generate unbalanced weight of the ascending/descending data in the inversion leading to possibly biased solutions [6]. We have therefore chosen to equally weight ascending and descending datasets by correcting covariance matrix from the datasets relative magnitude as expressed by Dumont et al. [6]: where C d i is the covariance sub-matrix of the i-th dataset, N d is the number of datasets, χ 2 0 i and χ 2 0 tot are the reference misfits (corresponding to a null displacement model) of the i-th dataset and of all datasets, respectively. In the Bayesian inference framework, this correction gives the same likelihood to each dataset.
The inversion procedures are conducted in 2 stages: a search and an appraisal stage. The search stage is done by combining forward model computations with a Neighborhood Algorithm [14], where the neighborhood is defined by Voronoi cells. Forward computations are conducted with a 3D Boundary element method which takes fractures, reservoirs and topographies into account. Source interactions are taken into account. Sources and topographies are meshed with triangular elements. To avoid edge effects, the topography mesh has an extension which is at least ten times larger than the intrusion [15]. Following previous estimations for Piton de la Fournaise [16,2], a Young's modulus of 5 GPa and a Poisson's coefficient of 0.25 were used. These values are consistent with in-situ measurements for basaltic volcanoes [17]. For each forward model, the displacements are computed for a unit over-pressure of 1 MPa which is then scaled to fit the data magnitude (see details in Tridon et al. [3]). If needed to improve the data fit, a shear stress component is added and scaled jointly with the over-pressure. Here the geometry used to represent sheet intrusion is approximated by a quadrangular dislocations defined by a set of parameters (see source descriptions in the next sections). We assumed the a priori probability density function to be uniform for each parameter.
The inversion algorithm starts by randomly picking N 1 points in the parameter space, where N 1 is an exponential function of the number of parameters, k, N 1 = 1.88 × e 0.65×k [3], which is determined by the number of natural neighbors [18]. For subsequent iteration, following Fukushima et al. [2], 50 points are drawn in the neighborhood of the 50 lowest misfit points. This enables the search to be explorative enough, while keeping the inversion time reasonable. Iterations stop when one of the following conditions is reached: (1) the standard deviation on χ 2 of the 50 th last model is smaller than 0.05, (2) the scaled model parameter values standard deviation of the 50 th last model is smaller than 0.05, (3) the maximum number of iterations have been reached (depending on the number of parameters considered).
In order to obtain numerically computable misfits, the interferograms are subsampled with a circular subsampling or a quadtree algorithm [19], to give a large number of points in high deforming areas, and a small number of points in low deforming areas.
In the appraisal stage of the inversions, 1D and 2D posterior probability functions are estimated using the Bayesian inference [20]. The model population calculated during the searching stage is resampled by a Monte Carlo integration, allowing posterior probability density functions to be reconstructed. The confidence interval and the mean model can then be derived from the 1D posterior probability density functions.

Goodness of fit estimation
We estimate the best model fit in a more comprehensive manner than misfit by defining the percentage of explained data: This estimator indicates how well the model fits the data relative to the data magnitude expressed by a null reference model (u T obs u obs ). It allows cross comparison of results between inversions having different number of datasets, data points or weighting. In addition, we compute the complementary RMS estimator giving an absolute mean value of the residuals of the models allowing to compare the residuals relatively to the measurement uncertainties.
where N is the number of data points.

Source type
Three different type of sources were used in this study. For intrusions that reached the surface with the fissure locations and orientations well defined, we used a quadrangular dislocation link to the surface by a given number of en echelon fractures set manually (Dyke source). For intrusions that failed to reach the surface or when fissures can not be defined, we used a quadrangular dislocation not connected to the ground surface (Deep Intrusion source). Finally, for intrusions that involves non-negligible shear component, we simplify the intrusion geometry by the same quadrangular dislocation as Deep Intrusion source but with an additional shear stress component which is linearly inverted as for the pressure value (Deep Sheared Intrusion source).
Dyke: quadrangular dislocation connected to the surface Deep sheared intrusion: quadrangular dislocation non-connected to the surface and submitted to shear stress Fig. S3: Geometric parameters used to define deep intrusion quadrangle geometries subjected to pressure and shear stress changes on a circular patch. Geometrical parameters are the same as deep intrusion quadrangle except the addition of the parameters Cx, Cy, r and rake which represent the x and y coordinates of the center of the circular patch, its radius and the rake of the shearing, respectively. Note that for all models done in this study, Cx, Cy, r are fixed to pressurized the whole the quadrangle. For more detail please refer to Tridon et al. [3] from which the figure is taken.  : Intrusive events since 1998 and available InSAR data. Each red patch indicates an intrusive event, with a width corresponding its duration. Events noted "d" or "a" indicate that only a descending or an ascending InSAR data is available, respectively, while "G" indicate that only GNSS data are available. The "x" notation indicates that a model has been determined for the eruption but not interpreted due to a too low signal to noise ratio. SAR data availability is given for each satellite: ENVISAT-ASAR (ENVISAT), Sentinel-1 (S1), RADARSAT-1 (RS1), RADARSAT-2 (RS2), COSMO-SkyMed (CSK), TerraSAR-X (TSX), TandDEM-X (TDX), ALOS-PALSAR-1 (ALOS1), ALOS-PALSAR-2 (ALOS2). Table S1: List of Piton de la Fournaise's eruptions and model's goodness of fit. Rift zone intruded at each event is indicated based on models. Volume (Vol), overpressure (P), shear stress (S), mean opening (Op) and shear (Sh) values of the intrusion are given, note that shear component is typically weak except for 4 events showing more than one meter of slip. Percentage of explained data (%ED, eq. 5) and RMS in meter (eq. 6) are given for ascending data (Asc), descending data (Dsc) and all data (tot